      subroutine sanm1n(il,ir)

      use param

#include <global.h>
      include 'ntbytes.h'
      include 'geom.h'
      include 'xsec.h'
      include 'sanm.h'
      include 'eigv.h'
      include 'sanmjin.h'

      integer,intent(in) :: il,ir
      real(NBF) cp011,cp012,cp021,cp022,cp111,cp112,cp121,cp122,cp211,cp212   &
    ,         cp221,cp222,cp311,cp312,cp321,cp322,cp411,cp412,cp421,cp422     &
    ,         krt11,krt12,krt21,krt22,rkrt11,rkrt12,rkrt21,rkrt22             &
    ,         ksq11,ksq12,ksq21,ksq22,rksq11,rksq12,rksq21,rksq22             &
    ,         skrt11,skrt12,skrt21,skrt22,ckrt11,ckrt12,ckrt21,ckrt22         &
    ,         eig11,eig12,eig21,eig22,beta11,beta12,beta21,beta22             &
    ,         rskrt11,rskrt12,rskrt21,rskrt22,cpt11,cpt12,cpt13,cpt14         &
    ,         cpt21,cpt22,cpt23,cpt24,cfm011,cfm012,cfm021,cfm022             &
    ,         cfm111,cfm112,cfm121,cfm122,csm11,csm12,csm21,csm22             &
    ,         error1,error2
!      real(NBF),dimension(:) :: j1_err(ng),j2_err(ng),j1_old(ng),j2_old(ng) 

!#define test
#ifdef test
      qs(0,1,1)=0
      qs(1,1,1)=0
      qs(2,1,1)=0
      qs(3,1,1)=0
      qs(4,1,1)=0
      qs(0,tnode,1)=0
      qs(1,tnode,1)=0
      qs(2,tnode,1)=0
      qs(3,tnode,1)=0
      qs(4,tnode,1)=0
      avgflx(1,1,1)=0.727729
      avgflx(1,tnode,1)=0.7277288
      avgflx(2,1,1)=0.04000650
      avgflx(2,tnode,1)=0.04000649      
      qs(0,1,2)=0.018880000
      qs(1,1,2)=0.034057800
      qs(2,1,2)=0.033608531
      qs(3,1,2)=0.046498862
      qs(4,1,2)=0.031526482
      qs(0,tnode,2)=0.018880000
      qs(1,tnode,2)=-0.034057800
      qs(2,tnode,2)=0.033608531
      qs(3,tnode,2)=-0.046498862
      qs(4,tnode,2)=0.031526482
      avgflx(1,1,2)=0.7039248
      avgflx(1,tnode,2)=0.7039246
      avgflx(2,1,2)=-0.001443321
      avgflx(2,tnode,2)=-0.001443322
#endif
      j11_old=0
      j12_old=0
      j21_old=0
      j22_old=0

      do iin=1,100
        fsc=0
!
        if (il.eq.1) then
          do im=1,ng
          i=1
! substitutions
          krt11=krt(1,i,im)
          krt21=krt(2,i,im)
          rkrt11=1/krt11
          rkrt21=1/krt21
          ksq11=ksq(1,i,im)
          ksq21=ksq(2,i,im)
          rksq11=rksq(1,i,im)
          rksq21=rksq(2,i,im)

          skrt11=sinh(krt11)
          skrt21=sinh(krt21)
          ckrt11=cosh(krt11)
          ckrt21=cosh(krt21)

          eig11=eigvec1(i,im)
          eig21=eigvec2(i,im)
          beta11=betax(1,i,im)
          beta21=betax(2,i,im)

          rskrt11=1/sinh(krt11)
          rskrt21=1/sinh(krt21)
! the coefficients of the particular solution 
! cp/order/moment/l-rnode/
          cp011=imq1(i,im)*(5*qs(4,i,im)*(21+2*ksq(1,i,im))+qs(0,i,im)*kqu(1,i,im) &
               +3*qs(2,i,im)*ksq(1,i,im))*rksx(1,i,im)
          cp021=imq2(i,im)*(5*qs(4,i,im)*(21+2*ksq(2,i,im))+qs(0,i,im)*kqu(2,i,im) &
               +3*qs(2,i,im)*ksq(2,i,im))*rksx(2,i,im)
          cp111=imq1(i,im)*(qs(1,i,im)*ksq(1,i,im)+15*qs(3,i,im))*rkqu(1,i,im)
          cp121=imq2(i,im)*(qs(1,i,im)*ksq(2,i,im)+15*qs(3,i,im))*rkqu(2,i,im)
          cp211=imq1(i,im)*(qs(2,i,im)*ksq(1,i,im)+35*qs(4,i,im))*rkqu(1,i,im)
          cp221=imq2(i,im)*(qs(2,i,im)*ksq(2,i,im)+35*qs(4,i,im))*rkqu(2,i,im)
          cp311=imq1(i,im)*qs(3,i,im)*rksq(1,i,im)
          cp321=imq2(i,im)*qs(3,i,im)*rksq(2,i,im)
          cp411=imq1(i,im)*qs(4,i,im)*rksq(1,i,im)
          cp421=imq2(i,im)*qs(4,i,im)*rksq(2,i,im)
! Coefficients of the Cosh term
          cosa(1,i,im)=-(rskrt11*krt11*(eig11*cp011-eig21*cp011 &
                     -avgflx(1,i,im)+eig21*avgflx(2,i,im)))/(eig11-eig21)      
         cosa(2,i,im)=(rskrt21*krt21*(eig21*cp021-avgflx(1,i,im) &
                     +eig11*(-cp021+avgflx(2,i,im))))/(eig11-eig21)
! Coefficients of the sinh term
          cfm011=-0.125*(5+4*eig11)*skrt11-ckrt11*krt11*beta11*(2+eig11)
          cfm012=-0.125*(5+4*eig21)*skrt21-ckrt21*krt21*beta11*(2+eig21)
          cfm021=0.125*skrt11*(eig11-5)-ckrt11*krt11*beta21
          cfm022=0.125*skrt21*(eig21-5)-ckrt21*krt21*beta21

          cpt11=cp111-3*cp211+6*cp311-10*cp411
          cpt12=cp121-3*cp221+6*cp321-10*cp421
          cpt13=-cp011+cp111-cp211+cp311-cp411
          cpt14=-cp021+cp121-cp221+cp321-cp421

          csm11=beta11*(cpt11*(eig11+2)+cpt12*(eig21+2))+0.625*(cpt13+cpt14) &
               +0.5*(eig11*cpt13+eig21*cpt14)-beta11*(cosa(1,i,im)*krt11*skrt11*(eig11+2) &
               +cosa(2,i,im)*krt21*skrt21*(eig21+2))-cosa(1,i,im)*ckrt11*(0.5*eig11+0.625) &
               -cosa(2,i,im)*ckrt21*(0.5*eig21+0.625) 
         csm12=0.125*(5*(cpt13+cpt14)-eig11*cpt13-eig21*cpt14+cosa(1,i,im)*ckrt11*(eig11-5) &
               +cosa(2,i,im)*ckrt21*(eig21-5))+beta21*(cpt11+cpt12-cosa(1,i,im)*krt11*skrt11 &
               -cosa(2,i,im)*krt21*skrt21)
          sinb(1,i,im)=(cfm022*csm11-cfm012*csm12)/(cfm011*cfm022-cfm012*cfm021)
          sinb(2,i,im)=(cfm021*csm11-cfm011*csm12)/(-cfm011*cfm022+cfm012*cfm021)
! Net current on the left boundary
          jnet(1,i,im)=-beta11*(krt11*(2+eig11)*(sinb(1,i,im)*ckrt11-cosa(1,i,im)*skrt11) &
                      +krt21*(2+eig21)*(sinb(2,i,im)*ckrt21-cosa(2,i,im)*skrt21)+eig11*cpt11 &
                      +eig21*cpt12+2*(cpt11+cpt12))
         jnet(2,i,im)=-beta21*(krt11*(sinb(1,i,im)*ckrt11-cosa(1,i,im)*skrt11) &
                      +krt21*(sinb(2,i,im)*ckrt21-cosa(2,i,im)*skrt21)+cpt11+cpt12)
! Surface flux
          sflux(1,i,im)=eig11*(ckrt11*cosa(1,i,im)-skrt11*sinb(1,i,im)-cpt13) &
                       +eig21*(ckrt21*cosa(2,1,im)-skrt21*sinb(2,i,im)-cpt14)
         sflux(2,i,im)=ckrt11*cosa(1,i,im)+ckrt21*cosa(2,i,im)-skrt11*sinb(1,i,im) &
                       -skrt21*sinb(2,i,im)-cpt13-cpt14

! 5 coefficients of the 4-th order flux
          cf(0,i,im)=cosa(1,i,im)*eig11*skrt11*rkrt11+cosa(2,i,im)*eig21*skrt21*rkrt21 &
                    +eig11*cp011+eig21*cp021 
         cf(1,i,im)=3*(sinb(1,i,im)*eig11*(ckrt11*rkrt11-skrt11*rksq11) &
                    +sinb(2,i,im)*eig21*(ckrt21*rkrt21-skrt21*rksq21)) &
                    +eig11*cp111+eig21*cp121
         cf(2,i,im)=5*(cosa(1,i,im)*eig11*rkrt11*(-3*ckrt11*rkrt11+3*skrt11*rksq11+skrt11) &
                    +cosa(2,i,im)*eig21*rkrt21*(-3*ckrt21*rkrt21+3*skrt21*rksq21+skrt21)) &
                    +eig11*cp211+eig21*cp221
         cf(3,i,im)=7*(sinb(1,i,im)*eig11*rksq11*(ckrt11*rkrt11*(15+ksq11)-15*skrt11*rksq11-6*skrt11) &
                    +sinb(2,i,im)*eig21*rksq21*(ckrt21*rkrt21*(15+ksq21)-15*skrt21*rksq21-6*skrt21)) &
                    +eig11*cp311+eig21*cp321
         cf(4,i,im)=9*(cosa(1,i,im)*eig11*rkrt11*(rksq11*(-5*ckrt11*rkrt11*(21+2*ksq11)+105*skrt11*rksq11 &
                    +45*skrt11)+skrt11)+cosa(2,i,im)*eig21*rkrt21*(rksq21*(-5*ckrt21*rkrt21*(21+2*ksq21) &
                    +105*skrt21*rksq21+45*skrt21)+skrt21))+eig11*cp411+eig21*cp421
! Update the 4-th order source coefficients
          do l=0,4
            fsc(l,i)=fsc(l,i)+xsnf(i,im)*cf(l,i,im)
          enddo
        enddo
!
        do im=1,ng
          scs=0  
          do m=1,ng   
            do l=0,4
              scs(l,i)=scs(l,i)+xssm(i,im)%from(m)*cf(l,i,m)      
            enddo
          enddo
          do l=0,4
            qs(l,i,im)=reigv*xchi(i,im)*fsc(l,i)+scs(l,i)          
          enddo
        enddo ! for im
!
        elseif (ir.eq.tnode) then
          do im=1,ng
! the coefficients of the particular solution 
! cp/order/moment/l-rnode/
          j=tnode
! substitutions
          krt12=krt(1,j,im)   
          krt22=krt(2,j,im)
          rkrt12=1/krt12
          rkrt22=1/krt22
          ksq12=ksq(1,j,im)
          ksq22=ksq(2,j,im)
          rksq12=rksq(1,j,im)
          rksq22=rksq(2,j,im)
   
          skrt12=sinh(krt12)      
          skrt22=sinh(krt22)      
          ckrt12=cosh(krt12)      
          ckrt22=cosh(krt22)      
                                  
          eig12=eigvec1(j,im) 
          eig22=eigvec2(j,im) 
          beta12=betax(1,j,im)
          beta22=betax(2,j,im)
                                  
          rskrt12=1/sinh(krt12)   
          rskrt22=1/sinh(krt22)   
!
          cp012=imq1(j,im)*(5*qs(4,j,im)*(21+2*ksq(1,j,im))+qs(0,j,im)*kqu(1,j,im) &
               +3*qs(2,j,im)*ksq(1,j,im))*rksx(1,j,im)
         cp022=imq2(j,im)*(5*qs(4,j,im)*(21+2*ksq(2,j,im))+qs(0,j,im)*kqu(2,j,im) &
               +3*qs(2,j,im)*ksq(2,j,im))*rksx(2,j,im)
          cp112=imq1(j,im)*(qs(1,j,im)*ksq(1,j,im)+15*qs(3,j,im))*rkqu(1,j,im)
          cp122=imq2(j,im)*(qs(1,j,im)*ksq(2,j,im)+15*qs(3,j,im))*rkqu(2,j,im)
          cp212=imq1(j,im)*(qs(2,j,im)*ksq(1,j,im)+35*qs(4,j,im))*rkqu(1,j,im)
          cp222=imq2(j,im)*(qs(2,j,im)*ksq(2,j,im)+35*qs(4,j,im))*rkqu(2,j,im)
          cp312=imq1(j,im)*qs(3,j,im)*rksq(1,j,im)
          cp322=imq2(j,im)*qs(3,j,im)*rksq(2,j,im)
          cp412=imq1(j,im)*qs(4,j,im)*rksq(1,j,im)
          cp422=imq2(j,im)*qs(4,j,im)*rksq(2,j,im)
! Coefficients of the cosh term
          cosa(1,j,im)=-(rskrt12*krt12*(eig12*cp012-eig22 &
                      *cp012-avgflx(1,j,im)+eig22*avgflx(2,j,im)))/(eig12-eig22)      
         cosa(2,j,im)=(rskrt22*krt22*(eig22*cp022-avgflx(1,j,im) &
                      +eig12*(-cp022+avgflx(2,j,im))))/(eig12-eig22)
! Coefficients of the sinh term
          cfm111=-0.125*(5+4*eig12)*skrt12-ckrt12*krt12*beta12*(2+eig12)
          cfm112=-0.125*(5+4*eig22)*skrt22-ckrt22*krt22*beta12*(2+eig22)
          cfm121=0.125*skrt12*(eig12-5)-ckrt12*krt12*beta22
          cfm122=0.125*skrt22*(eig22-5)-ckrt22*krt22*beta22

          cpt21=cp112+3*cp212+6*cp312+10*cp412
          cpt22=cp122+3*cp222+6*cp322+10*cp422
          cpt23=cp012+cp112+cp212+cp312+cp412
          cpt24=cp022+cp122+cp222+cp322+cp422

          csm21=beta12*(cpt21*(eig12+2)+cpt22*(eig22+2))+0.625*(cpt23+cpt24) &
               +0.5*(eig12*cpt23+eig22*cpt24)+beta12*(cosa(1,j,im)*krt12*skrt12*(eig12+2) &
               +cosa(2,j,im)*krt22*skrt22*(eig22+2))+cosa(1,j,im)*ckrt12*(0.5*eig12+0.625) &
               +cosa(2,j,im)*ckrt22*(0.5*eig22+0.625)
         csm22=-0.125*(eig12*cpt23+eig22*cpt24-5*(cpt23+cpt24)+cosa(1,j,im)*ckrt12*(eig12-5) &
               +cosa(2,j,im)*ckrt22*(eig22-5))+beta22*(cpt21+cpt22+cosa(1,j,im)*krt12*skrt12 &
               +cosa(2,j,im)*krt22*skrt22)
        
          sinb(1,j,im)=(cfm122*csm21-cfm112*csm22)/(cfm111*cfm122-cfm112*cfm121)
          sinb(2,j,im)=(cfm121*csm21-cfm111*csm22)/(-cfm111*cfm122+cfm112*cfm121)
! Net current on the right boundary
          jnet(1,j+1,im)=-beta12*(krt12*(2+eig12)*(sinb(1,j,im)*ckrt12+cosa(1,j,im)*skrt12) &
                        +krt22*(2+eig22)*(sinb(2,j,im)*ckrt22+cosa(2,j,im)*skrt22) &
                        +eig12*cpt21+eig22*cpt22+2*(cpt21+cpt22))
         jnet(2,j+1,im)=-beta22*(krt12*(sinb(1,j,im)*ckrt12+cosa(1,j,im)*skrt12) &
                        +krt22*(sinb(2,j,im)*ckrt22+cosa(2,j,im)*skrt22)+cpt21+cpt22)
! Surface flux
          sflux(1,j+1,im)=eig12*(ckrt12*cosa(1,j,im)+skrt12*sinb(1,j,im)+cpt23) &
                         +eig22*(ckrt22*cosa(2,j,im)+skrt22*sinb(2,j,im)+cpt24)
         sflux(2,j+1,im)=ckrt12*cosa(1,j,im)+ckrt22*cosa(2,j,im)+skrt12*sinb(1,j,im) &
                         +skrt22*sinb(2,j,im)+cpt23+cpt24

! 5 coefficients of the 4-th order flux
          cf(0,j,im)=cosa(1,j,im)*eig12*skrt12*rkrt12+cosa(2,j,im)*eig22*skrt22*rkrt22 &
                    +eig12*cp012+eig22*cp022
         cf(1,j,im)=3*(sinb(1,j,im)*eig12*(ckrt12*rkrt12-skrt12*rksq12) &
                    +sinb(2,j,im)*eig22*(ckrt22*rkrt22-skrt22*rksq22)) &
                    +eig12*cp112+eig22*cp122
         cf(2,j,im)=5*(cosa(1,j,im)*eig12*rkrt12*(-3*ckrt12*rkrt12+3*skrt12*rksq12+skrt12) &
                    +cosa(2,j,im)*eig22*rkrt22*(-3*ckrt22*rkrt22+3*skrt22*rksq22+skrt22)) &
                    +eig12*cp212+eig22*cp222
         cf(3,j,im)=7*(sinb(1,j,im)*eig12*rksq12*(ckrt12*rkrt12*(15+ksq12)-15*skrt12*rksq12-6*skrt12) &
                    +sinb(2,j,im)*eig22*rksq22*(ckrt22*rkrt22*(15+ksq22)-15*skrt22*rksq22-6*skrt22)) &
                    +eig12*cp312+eig22*cp322
         cf(4,j,im)=9*(cosa(1,j,im)*eig12*rkrt12*(rksq12*(-5*ckrt12*rkrt12*(21+2*ksq12)+105*skrt12*rksq12 &
                    +45*skrt12)+skrt12)+cosa(2,j,im)*eig22*rkrt22*(rksq22*(-5*ckrt22*rkrt22*(21+2*ksq22) &
                    +105*skrt22*rksq22+45*skrt22)+skrt22))+eig12*cp412+eig22*cp422
! Update the 4-th order source coefficients
          do l=0,4
            fsc(l,j)=fsc(l,j)+xsnf(j,im)*cf(l,j,im)
          enddo
        enddo ! for im
!
          do im=1,ng
            scs=0  
            do m=1,ng   
              do l=0,4
                scs(l,j)=scs(l,j)+xssm(j,im)%from(m)*cf(l,j,m)      
              enddo
            enddo  
            do l=0,4
              qs(l,j,im)=reigv*xchi(j,im)*fsc(l,j)+scs(l,j)             
            enddo
          enddo
        endif
! Current Convergence condition
        if (bcb.eq.1.and.il.eq.1) then
          k=i
        elseif (bcu.eq.1.and.ir.eq.tnode) then
          k=j+1
        endif

        do im=1,ng
          j1_err(im)=abs(jnet(1,k,im)-j1_old(im))
          j2_err(im)=abs(jnet(2,k,im)-j2_old(im))
          j1_old(im)=jnet(1,k,im)
          j2_old(im)=jnet(2,k,im)
        enddo
        error1=maxval(j1_err)
        error2=maxval(j2_err)
        if (error1.lt.epsm6.and.error2.lt.epsm6) then
         print*,"node=",k,"nodal=",iin
         exit
        endif 
      enddo   ! for iin

      end subroutine
